Genome-wide chromatin contacts of super-enhancer-associated lncRNA identify LINC01013 as a regulator of fibrosis in the aortic valve

Calcific aortic valve disease (CAVD) is characterized by a fibrocalcific process. The regulatory mechanisms that drive the fibrotic response in the aortic valve (AV) are poorly understood. Long noncoding RNAs derived from super-enhancers (lncRNA-SE) control gene expression and cell fate. Herein, multidimensional profiling including chromatin immunoprecipitation and sequencing, transposase-accessible chromatin sequencing, genome-wide 3D chromatin contacts of enhancer-promoter identified LINC01013 as an overexpressed lncRNA-SE during CAVD. LINC01013 is within a loop anchor, which has contact with the promoter of CCN2 (CTGF) located at ~180 kb upstream. Investigation showed that LINC01013 acts as a decoy factor for the negative transcription elongation factor E (NELF-E), whereby it controls the expression of CCN2. LINC01013-CCN2 is part of a transforming growth factor beta 1 (TGFB1) network and exerts a control over fibrogenesis. These findings illustrate a novel mechanism whereby a dysregulated lncRNA-SE controls, through a looping process, the expression of CCN2 and fibrogenesis of the AV.


Introduction
Calcific aortic valve disease (CAVD) is the most frequent heart valve disorder [1]. The progression of CAVD culminates in aortic valve (AV) stenosis, which can be treated by aortic valve replacement (AVR) in symptomatic patients. AVR either surgical or by percutaneous approach is associated with elevated cost and significant morbidity/mortality [2]. Hence, the identification of key underpinning processes that control the development of CAVD could lead to medical therapy to prevent the progression of CAVD. Analyses of surgically explanted AVs have shown that the thickening of leaflets, a hallmark feature of CAVD, is associated with an excess of extracellular matrix components such as collagen fibrils [2]. Research conducted over the last several years has underlined several molecular processes involved in the development of CAVD and has highlighted that valve interstitial cells (VICs) have a high plasticity [3]. As such, according to different cues VICs may transition from quiescent to activated cells with an increased potential for fibrogenesis [3]. Studies have consistently shown that the transforming growth factor (TGF) beta pathway was involved in fibrogenesis and the development of CAVD [4]. Pre-clinical and clinical data have highlighted the benefit of targeting the TGF beta pathway to reduce the cardiac fibrosis [5] and the key role of the TGF beta-CCN2/CTGF axis [6]. However, key regulators of the TGF beta pathway that orchestrate and instruct a fibrogenic programme in the AV are still largely unknown.
Transcription factors and enhancers control gene expression pattern [7]. Super-enhancers (SE) are groups of distant-acting cis-regulatory elements (enhancers) in close proximity, which are enriched in chromatin loops with the promoter of genes [8]. The folding pattern of the chromatin brings distant elements, such as SE, in contact with gene promoters and exert a control on the transcriptional process. Studies have highlighted that SE are cell-specific and control cell fate and identity [8]. Enhancer and SE may express noncoding RNA including long noncoding RNA (lncRNA) [9,10]. The expression of lncRNA derived from SE (lncRNA-SE) provides another layer of control on gene expression [10]. Studies on cardiac fibrosis have demonstrated that lncRNA or lncRNA-SE could be targeted to reduce cardiac fibrosis [11,12]. In this work, we hypothesized that dysregulated lncRNA-SE expressed by VICs may contribute to fibrogenesis, an early and key process in the development of CAVD. By using functional assays and genomic multidimensional profiling of VICs and AVs (S1 Fig), we identified LINC01013, a CAVD-dysregulated lncRNA-SE within a regulatory DNA loop that controls the transcriptional process of gene enriched in the TGF beta pathway and promotes fibrogenesis.

Super-enhancer in valve interstitial cells
SEs exhibit a strong enrichment for the active H3K27ac regulatory mark and play a major role in the signaling pathways. We thus wondered if SEs could be involved in the biology of the AV. For that purpose, in VICs, we first generated chromatin immunoprecipitation for H3K27ac followed by sequencing to identify active regulatory elements. We identified 82,662 strong peaks corresponding to putative regulatory elements. We next used these data to identify SEs by using the ranking of super-enhancer (ROSE) algorithm as described previously [13,14]. ROSE is the most utilized approach to identify SEs; briefly the algorithm performs a ranking of strong and close H3K27ac peaks, and then identifies SEs from peaks plotted on a scatter plot upper an inflection point for which the slope is 1. Thus, in VICs, we identified 1085 SEs (S1 Table). We wanted to replicate these data by using another algorithm integrated into the HOMER tools, and we observed that 83.3% of the SEs identified by HOMER (S2 Table) are also found by the ROSE algorithm. By using GREAT (see the Methods section), we found that the 1085 SEs were highly enriched in gene ontology (GO) for extracellular matrix organization (GO:0030198) (P = 2.01 x 10 −21 ) (Fig 1A). Active regulatory elements like enhancers and SEs are characterized by a chromatin accessibility. In VICs, we performed an assay for transposase-accessible chromatin and massively parallel sequencing (ATAC-seq) to capture the open chromatin regions. We found that there was a significant and a positive correlation of the ATAC-seq signal close to the SEs (observed/expected ratio = 5.2, P<1 x 10 −16 binomial test) ( Fig 1B): these findings attest the importance of the SEs identified as active regulatory elements. By using an extensive annotation of lncRNAs (GENCODE [15] version 32) that was intersected with the SE dataset, we identified 324 SEs from which at least one lncRNA is transcripted, called lncRNA-SE (S3 Table), which represent 30% of all the SE in VICs ( Fig  1C). Based on the ROSE algorithm [8], lncRNA-SE were ranked higher than regular SE (R-SE), from which no lncRNA is transcribed (P<0.0001, Wilcoxon ranked-sum test) ( Fig  1D). We thus wondered whether lncRNA-SE were associated with increased open chromatin compared to R-SE. In ATAC-seq, the tag density for the core region (±5kB) was significantly increased in lncRNA-SE compared to R-SE (P<0.0001, Wilcoxon ranked-sum test) (Fig 1E  and 1F). SE are enriched in 3D chromatin looping involved in gene regulation [14]. In VICs, genome-wide enhancer-promoter 3D mapping was evaluated by using H3K27ac HiChIP, which generated 454,964,884 unique paired-end tags. Fig 2A shows an example of an interaction matrix for the chromosome 6 with a resolution of up to 1kb. S2 Fig shows interaction matrices for the others chromosome. In total, in VICs we mapped 36,229 high-confidence loops at false discovery rate (FDR)<0.01. The mean number of protein coding gene promoters (±1 kb from transcription start site) mapped by 3D chromatin looping was 1.2 and 1.9 for R-SE and lncRNA-SE respectively. Taken together, these data highlight that VIC-associated lncRNA-SE are highly ranked SE with increased open chromatin and are spatially coordinated with distant protein coding genes.

Differentially expressed lncRNA-SE and 3D genome mapping
We next wondered if lncRNA-SE were differentially expressed during the development of CAVD. By using RNA-sequencing (RNA-seq) in 10 control nonmineralized AVs and 19 surgically-explanted mineralized AVs (CAVD), we profiled the expression of lncRNAs. RNA sequencing was annotated with GENCODE version 32. Fig 2B presents a heatmap of differentially expressed lncRNAs between control nonmineralized and mineralized AVs (CAVD). When comparing CAVD to control AVs, 178 and 111 lncRNAs were upregulated and downregulated respectively (FDR<0.05) (S4 Table). Among the differentially regulated lncRNAs there were 18 lncRNA-SE (10 upregulated and 8 downregulated in CAVD) (S5 Table). In order to identify the putative genes regulated by lncRNA-SE we mapped these regulatory elements to their protein coding genes (±1 kb from the transcription start site) by using chromatin looping information obtained with the H3K27ac HiChIP in VICs. Among the 18

Fig 1. Super-enhancer (SE) in valve interstitial cells (VICs).
A) Gene ontology enrichment for SE in VICs (adjusted P-value). B) ATAC-seq signal relative to distance from SE. C) Pie chart showing that 30% of VICs SE intersect with lncRNA (lncRNA-SE); SEs for which no lncRNA is transcribed are defined as regular SE (R-SE). D) The ranking of super-enhancer (ROSE) algorithm indicates that lncRNA-SE rank higher than R-SE in VICs. E) Signal intensity for open chromatin (ATAC-seq) was higher at the core of lncRNA-SE compared with R-SE. F) Box plot revealing higher ATAC-seq signal for lncRNA-SE. https://doi.org/10.1371/journal.pgen.1010010.g001

PLOS GENETICS
Super-enhancer-associated lncRNA and fibrogenesis of the aortic valve

PLOS GENETICS
Super-enhancer-associated lncRNA and fibrogenesis of the aortic valve differentially expressed lncRNA-SE, we found that in total 4 were having significant chromatin contacts with distant promoters of protein coding genes (S6 Table). The mean distance between a differentially expressed lncRNA-SE and a protein coding gene was~89 kb. We wondered whether genes mapped by H3K27ac HiChIP and having loops with lncRNA-SE could be part of the GO for extracellular matrix organization, the highest enriched ontology for SE in VICs. Fig 2C shows a Venn diagram of genes mapped by chromatin contact with differentially expressed lncRNA-SE and genes encompassed in the GO for extracellular matrix organization. One gene promoter for CCN2 (6q23.2) having loops with a differentially expressed lncRNA-SE overlapped with the GO for extracellular matrix organization. We next examined the genome 3D interactions at 6q23.2 by using virtual 4C (v4C), which provides an anchor as a viewpoint of all interactions (with the anchor) within a genomic region and visualized in 2D. Fig 2D  shows v4C representation of chromatin interactions at 6q23.2 in VICs where LINC01013/ LOC100507254, a poorly conserved lncRNA-SE according to the NONCODE database [16], interacts with the promoter region of CCN2 located~180 kb upstream. High confidence significant loops identified strong chromatin contacts between the promoter of LINC01013 and the promoter of CCN2 as well as many significant non-coding chromatin contacts at these loci. LINC01013 is overexpressed in mineralized compared to control nonmineralized AVs and SE region for LINC01013 extends over~30 kb. CCN2 encodes for the connective tissue growth factor (CTGF), a matricellular protein involved in cell adhesion, proliferation and synthesis of extracellular matrix [17,18]. Interrogation of enhancer-promoter chromatin contacts using H3K27ac HiChIP at the 6q23.2 locus showed that similarly to VICs, LINC01013 interacts with CCN2 in human coronary artery smooth muscle cells (HCASMC), whereas there is no chromatin looping at this locus in blood cells such as naive T cells and GM12878 (B cell line) ( Fig 2D). These data underline that spatial organization of chromatin at the LINC01013 locus is likely specific to mesenchymally-derived cells. The S3 Fig gives a summary of the sequential analysis which conducted to the prioritization of LINC01013.

Characterization of the LINC01013 locus
We confirmed with a higher resolution the data obtained in HiChIP by using chromatin conformation capture (3C) and pairs of primers covering the region of LINC01013 (S7 Table). In VICs, the 3C assay showed that the promoter region of CCN2 (used as an anchor) interacted with the promoter region of LINC01013 ( Fig 3A). Alternative splicing of LINC01013 produces transcripts with 2 or 3 exons. In VICs, only the 3 exons transcript was expressed, whereas the 2 exons transcript was not detected by RT-qPCR ( Fig 3B). We generated data to confirm the RNA-seq analysis. We measured the level of transcripts encoding for LINC01013 in 18 control nonmineralized and 26 mineralized AVs by using quantitative RT-qPCR. Table 1 presents the demographic data for the patients for whom the AVs were obtained. We found that LINC01013 was increased by 2.95-fold in mineralized compared to control AVs (P<0.0001; age-and sex-adjusted P = 0.0003) (Fig 3C). The noncoding potential of LINC01013 was confirmed by using CPAT [19] (coding probability = 0.04), a logistic regression algorithm to probe the coding potential (S8 Table). In VICs, a ChIP-seq with an antibody directed against H3K4me1, a histone mark of enhancer, showed strong peaks over the LINC01013 region, which overlapped with H3K27ac and open chromatin ( Fig 2D). We next asked whether the dysregulation of LINC01013 in CAVD could be linked to a disease-associated epigenetic process. We measured by quantitative ChIP the enrichment of H3K27ac at the promoter of LINC01013 in control and mineralized AVs. Compared to control AVs, we found that the H3K27ac mark at the promoter of LINC01013 was increased by 5.41-fold in mineralized AVs ( Fig 3D). Taken together, these orthogonal data suggest that LINC01013 is located in extended histone marks of active regulatory elements and is spatially coordinated with CCN2.

LINC01013 and transcriptional response at CCN2
We assessed the expression of LINC01013 in the cytosol and the nuclear compartments of VICs. After cell fractionation, we observed that LINC01013 was expressed in both the nuclear

PLOS GENETICS
and cytosol fractions ( Fig 4A). To explore the functional role, we designed an antisense oligonucleotide (ASO) targeting LINC01013. In VICs, the transfection of ASO provided a strong reduction of LINC01013 ( Fig 4B). Transfection of the ASO against LINC01013 in VICs also led to a significant reduction of transcripts encoding for CCN2 measured after 72 hours ( Fig 4C). In the same line, the knockdown of LINC01013 was followed by a significant reduction of CCN2 (CTGF) at the protein level ( Fig 4D). Promoter-proximal pausing of RNA polymerase II (RNAPII) is a widespread regulatory mechanism of the transcriptional process [20]. NELF (negative elongation factor) is a protein complex interacting with the RNAPII [21] which restrains its progression usually 20 to 60 bp downstream from the transcription start site (TSS) [20]. The recruitment of the NELF complex is associated with the phosphorylation on serine 5 of the c-terminal domain of RNAPII (RNAPII-Ser5p) enriched at paused promoters [20,22]. Previous studies have identified that some noncoding RNAs can modulate the recruitment of the NELF complex by sequestering the NELF-E subunit [23,24]. Considering that LINC01013 has a nuclear localization and its locus is spatially coordinated with the promoter of CCN2, we hypothesized that LINC01013 can act as a decoy factor for the NELF-E subunit and thus modulate its recruitment at CCN2. We performed a RNA immunoprecipitation (RIP) assay to evaluate the possible interaction between LINC01013 and NELF-E. In VICs, a pull-down of the NELF-E subunit followed by a quantitative RT-qPCR showed that LINC01013 was enriched by 6-fold compared to the IgG control ( Fig 4E). In addition, quantitative ChIP demonstrated that the knockdown of LINC01013 in VICs led to a significant enrichment of the NELF-E subunit at the TSS of CCN2 ( Fig 4F). However, at the promoter of LINC01013, we observed a poor enrichment of NELF-E compared to the IgG control and no impact on the recruitment after the knockdown (S4A Fig). In line with these findings, we found that the knockdown of LINC01013 was accompanied by a higher level of RNAPII-Ser5p at the TSS of CCN2 ( Fig 4G). However, total RNAPII revealed a similar occupancy at the TSS of CCN2 with or without the knockdown ( Fig 4H). These data thus indicate that LINC01013 controls the magnitude of the

PLOS GENETICS
Super-enhancer-associated lncRNA and fibrogenesis of the aortic valve

PLOS GENETICS
Super-enhancer-associated lncRNA and fibrogenesis of the aortic valve transcriptional response at CCN2 through a modulation of RNAPPII pausing by acting as a decoy factor for the NELF-E subunit.

TGF beta pathway and LINC01013
We extracted a protein-protein interaction (PPI) network from RNA-seq in calcified AVs in order to evaluate the interaction of CCN2 and its pathway. For that purpose, we used an experimentally validated dataset including 156,317 PPIs. The extracted network in CAVD included 446 nodes (gene derived proteins) and 591 edges (connections) and was enriched for TGF beta signaling pathway (P adjusted = 4.17 x 10 −72 ) (Fig 5A and 5B). Highly connected nodes such as SMAD2, SMAD3, TGFB1 and COL1A1 were linked to CCN2. We measured in VICs the expression of COL1A1, a widely expressed collagen in the AV, in response to a knockdown of LINC01013 knowing its impact on the expression of CCN2. In VICs, the knockdown of LINC01013 lowered the protein level of COL1A1 by 55% after 72 hours ( Fig 6A). Also, the transfection of an ASO targeting LINC01013 in VIC cultures lowered the extracellular secretion of pro-collagen 1-α by 31% ( Fig 6B). Immunostaining of COL1A1 in VICs in culture demonstrated similar results on whole cells (S4B Fig). TGF beta is an enriched pathway in CAVD and an important promoter of VIC activation [3]. We tested whether TGFB1 may impact on the expression of CCN2 through the activation of LINC01013. First, in TGFB1-treated VIC cultures (10 ng/ml), we found that the expression of LINC01013 was increased by 1.9-fold ( Fig 6C). Similarly, the expression of CCN2 was increased by 3.1-fold in TGFB1-treated VIC cultures ( Fig 6D). Of interest, the ASO against LINC01013 abrogated the TGFB1-induced rise of CCN2 ( Fig 6D). Also, in western blot, the level of CCN2 increased following a treatment with TGFB1, whereas the knockdown of LINC01013 largely abrogated this rise ( Fig 6E). We next treated VIC cultures with TGFB1 and we measured the level of collagen with and without a knockdown of LINC01013. The knockdown of LINC01013 in VIC cultures prevented the TGFB1-induced rise of COL1A1 transcript and the secretion of pro-collagen type 1-α (Fig 6F and 6G). Data obtained from an immunostaining of COL1A1 in whole cells are consistent with these results (S4 Fig). Similarly, the knockdown of CCN2 by an ASO, which significantly reduced the transcript level of the target, also abrogated TGFB1-induced rise of COL1A1 (Fig 6H and 6I). Taken together, these data obtained in VICs indicate that TGFB1-mediated expression of collagen relies on a LINC01013-CCN2 pathway.

Discussion
Herein, genomic profiling of VICs identified that SE are enriched for extracellular matrix organization. Among these clusters of regulatory elements, lncRNA-SE are enriched in open chromatin and loops with gene promoters. We identified a lncRNA-SE, LINC01013, which is upregulated during CAVD. Functional assessment showed that LINC01013 is within a DNA loop that has contacts with the promoter of CCN2. LINC01013 is acting as a decoy factor for NELF-E, whereby it controls the transcription at CCN2. Increased expression of LINC01013 promotes the dissociation of NELF-E from chromatin and the expression of CCN2, which encodes for a matricellular protein that stimulates fibrogenesis (Fig 6J).

LncRNA and super-enhancers in VICs
In CAVD, previous studies have identified dysregulated lncRNAs such as TUG1 [25] and MALAT1 [26], which act as sponges of miR-204, a microRNA that regulates the level of the osteogenic transcription factor RUNX2. Also, the expression of the lncRNA H19 is increased in VICs during CAVD and controls the expression of NOTCH1, a key regulator of cell fate [27]. In the heart, CARMEN is a lncRNA-SE, which regulates cardiac specification of precursor cells [28]. To our knowledge, the present work is the first to identify and characterize dysregulated lncRNA-SE in CAVD. In VICs, mapping of SE revealed an enrichment for extracellular matrix organization, cell substrate adhesion and actin cytoskeleton. These functional enrichments represent salient biologic features of the AV, which are dysregulated during the These findings suggested that lncRNA-SE may participate to cell fate by controlling gene expression in regulatory DNA loops. LINC01013 was prioritized for further functional assessment as this lncRNA-SE has chromatin contacts with the promoter of CCN2, a gene involved in extracellular matrix organization. In explanted AVs, we found in CAVD that the promoter of LINC01013 was enriched for H3K27ac, a mark of active chromatin. Thus, these findings highlighted that lncRNA-SE are dysregulated during CAVD and are involved into several processes pertaining to matrix organization and activation of VICs.

LINC01013 regulates the transcription of CCN2
Transcription is a highly organized process involving the recruitment of regulatory complexes along with the RNAPII [29]. Distant-acting elements modulate the transcription process through DNA loops [30]. Herein, we found that LINC01013 is a decoy factor for NELF-E and exerts, through a long-distance looping, a control on the expression of CCN2. The NELF complex is known to be required for the phosphorylation on serine 5 of RNAPPII associated with RNAPII pausing [20,22]. Consistently, following the knockdown of LINC01013, we found that the level of RNAPII-Ser5p at the TSS of CCN2 was increased. Hence, during CAVD, higher expression of LINC01013 promotes the dissociation of NELF-E from chromatin and increases the transcriptional response of CCN2. Interestingly, the importance of LINC01013 in the transcriptional control of CCN2 is also observed in the endothelium. In human umbilical endothelial cells (HUVEC), the presence of SEs at the LINC01013 locus has been also reported as well as a strong interaction between the LINC01013 and the CCN2 locus [31]. Also, single-cell RNA sequencing experiments in HUVEC revealed that a sub-population of cells co-expressed LINC01013 and CCN2 [32]. In this regard, it appears that an unexplored role for LINC01013 on the expression of CCN2 in endothelial cells of the aortic valve might be an additional process involved into the pathobiology of CAVD.

TGF beta pathway
Transition of VICs towards activated cells is promoted by TGFB1. In a PPI network, CCN2 has several connections with different members of the TGF beta pathway. In line with previous works, we found in VICs that TGBFB1 induced the expression and secretion of collagen type 1, a predominant collagen within the AV. Of interest, we found that TGFB1 also induced the expression of LINC01013. We thus hypothesized that LINC01013 may regulate TGFB1-dependent fibrogenesis. To this effect, ASO-mediated knockdown of LINC01013 prevented TGFB1-induced expression of COL1A1 and secretion of pro-collagen type 1-α. Also, TGFB1-induced expression of COL1A1 was abrogated by the silencing of CCN2. Taken together, the present findings in VICs suggest that TGFB1-induced synthesis of collagen relied on a LINC01013-CCN2 pathway. These data are consistent with the role of CCN2 in fibrotic diseases.

Clinical implications
The present work underlined that LINC01013 regulates the transcription of CCN2, which promotes fibrogenesis, an underpinning process in the development of CAVD. Hence, these findings provide evidence that LINC01013 could be targeted in order to control CAVD. Further work on the role of key regulators of fibrogenesis in the AV could lead to novel therapies.

Limitations
Though we showed that LINC01013 is a transcriptional regulator of CCN2 by acting as decoy factor for NELF-E, the present study did not investigate the possible direct interaction of LINC01013 with DNA at this locus. Further experiments are needed to address the interaction of LINC01013 with chromatin. The expression of lncRNAs was documented in surgically explanted AVs and thus it was not possible to assess whether LINC01013 is expressed early during the development of CAVD. The role of LINC01013 was examined in vitro and whether targeting its expression in vivo would result in lower fibrosis of the AV is presently unknown, and the development of a mouse model for translational studies could constitute a challenge: the poor conservation of LINC01013 through species requires more investigation to identify a putative orthologue in mouse. However, conservation of synteny between human and mouse at the CCN2 and LINC01013 loci would suggest the presence of an orthologue, for which the sequence remains to be determined by using unbiased methods like rapid amplification of cDNA ends followed by long-read RNA sequencing (RACE-seq).

Conclusions
We provide evidence that lnRNA-SE are dysregulated during CAVD and participate to cell fate and extracellular matrix organization. Specifically, we identified LINC01013, a lnRNA-SE that is spatially coordinated with CCN2 and regulates the expression of CCN2 and the fibrotic response. Increased expression of LINC01013 promotes the expression of CCN2 and the synthesis of collagen. Targeting LINC01013 could prevent fibrosis of the AV.

Ethics statement
The protocol was approved by the ethical committee of the Institut Universitaire de Cardiologie et de Pneumologie de Québec. Written formal consent was obtained from the subjects.

VICs isolation
VICs were isolated from non-mineralized aortic valves obtained from patients during heart transplantation procedure. Aortic leaflets cut into pieces were incubated in 0.3% type I collagenase (Invitrogen, Thermo Fisher Scientific, ON, Canada) at 37˚C for 45 minutes and then filtered through a 70 μm mesh before a centrifugation for 5 minutes at 1,500 rpm. Cells were resuspended in complete media (DMEM, 10% FBS with L-glutamine and sodium pyruvate) and used between passages 3 to 7.

ChIP-seq library preparation and sequencing
VICs isolated were fixed with 1% formaldehyde at room temperature for 10 minutes and then quenched with glycine to 0.125 M. Cells were harvested in 1 mL of PBS 1X/5% BSA (Sigma-Aldrich, ON, Canada) and then centrifuged 5 minutes at 800 x g at 4˚C before removing the supernatant. Cells were resuspended in 1 mL of lysis buffer L1A/B (10 mM Hepes-KOH pH 7.5, 85 mM KCl, 1 mM EDTA, 0.5% NP-40 containing PMSF, sodium butyrate and PIC (Sigma-Aldrich, ON, Canada)) and incubated on ice for 10 minutes. Cellular extracts were centrifuged for 5 minutes at 8,000 x g, then nuclei were sonicated on ice with a Bioruptor Standard Sonicator (Diagenode, NJ, USA) in 450 μL of L3 buffer (5 mM Tris-HCl pH 7.5, 100 mM NaCl, 10 mM EDTA, 0.5 mM EGTA, 0.1% deoxycholate, 0.5% sarkosyl, supplemented with PMSF, sodium butyrate and PIC) for 3 minutes cycles alternating 10 seconds ON and 10 seconds OFF with 30% amplitude. Triton-X100 to final concentration of 1% was added and nuclear extracts were centrifugated for 5 minutes at 18

ATAC-seq library preparation and sequencing
Nuclei from approximatively 2.5 x 10 5 to 4 x 10 5 of VICs were treated as previously [33] described with some modifications. Transposition reactions were performed using the Tn5 Transposase and TD reaction buffer from the Nextera DNA library preparation kit (

ChIP-seq and ATAC-seq analysis
FASTQ files were mapped to the UCSC genome build hg19 using Bowtie2 with default settings. BigWIG and peak calling were respectively generated with bamcoverage tool and MACS2 with default parameters.

Super-enhancers analysis
Super-enhancers identification was performed as described previously [13,14] with a few modifications. First, MACS2 was used to identify enhancers based on H3K27ac peaks with a qvalue threshold of 1 x 10 −5 . H3K27ac peaks which overlapped with the ENCODE [34] blacklisted genomic regions were removed. The remaining H3K27ac peaks defined as enhancers were then used to identify super-enhancers (SE) by performing the ROSE [8] (Rank Ordering of Super-Enhancers) algorithm with default settings (stitching distance of 12.5 kb without TSS exclusion). GREAT [35] was used to determine Gene Ontology (GO) enrichment of SE; GREAT uses a binomial test to identify significant enrichments from the deviations of a theoretically expected distribution of observations to control for false positives, and gives a binomial p-value associated with the enrichments. SE ranking by ROSE were compared with a Wilcoxon ranked-sum test. Open chromatin enrichment for SE was evaluated with the Geno-metricCorr [36] package. lncRNA-SE and R-SE were defined from SE by using the Bedtools intercept function with the TSS of lncRNAs from the GENCODE [15] version 32 dataset. ATAC-seq signal around lncRNA-SE and R-SE were calculated with HOMER [37] and plot with R by using a loess smoothing method. The 5kb centered ATAC-seq signal from lncRNA-SE and R-SE was compared with a Wilcoxon ranked-sum test.

H3K27ac HiChIP
H3K27ac HiChIP data for GM12878, Naïve T cells and HCAMSC were used from GEO (accession number GSE101498) [38]. H3K27ac HiChIP for VICs was performed as described by Mumbach et al. [39]. About  , and 660 μL of water. Nuclei were then isolated by centrifugation at 2,500 x g for 5 minutes and resuspended in Nuclear lysis buffer (50 mM Tris-HCl pH 7.5, 10 mM EDTA, 1% SDS, supplemented with PIC, Sigma-Aldrich, ON, Canada) before sonication on ice with a Bioruptor Standard Sonicator (Diagenode, NJ, USA). Sample was clarify for 15 minutes at 16,100 x g at 4˚C and 800 μL of ChIP Dilution buffer (0.01% SDS, 1.1% Triton X-100, 1.2 mM EDTA, 16.7 mM Tris-HCl pH 7.5, 167 mM NaCl) were added. 34 μL of protein G Dynabeads (Life Technologies, Thermo Fisher Scientific, ON, Canada) were washed and resuspended in 50 μL of ChIP Dilution buffer. Sample was then added to beads, followed by an incubation at 4˚C for 1 hour with rotation. Sample was put on a magnet and supernatant removed into a new tube. 4 μg of H3K27ac antibody (#ab4729, Abcam, ON, Canada) were added, and the sample was incubated overnight at 4˚C with rotation. 34 μL of protein G Dynabeads in 50 μL of ChIP dilution buffer were added to the sample for an incubation at 4˚C for 2 hours. Beads were then washed sequentially with Low Salt Wash buffer (0.1% SDS, 1% Triton X-100, 2 mM EDTA, 20 mM Tris-HCl pH 7.5, 150 mM NaCl), High Salt Wash buffer (0.1% SDS, 1% Triton X-100, 2 mM EDTA, 20 mM Tris-HCl pH 7.5, 500 mM NaCl) and LiCl wash buffer (10 mM Tris-HCl pH 7.5, 250 mM LiCl, 1% NP-40, 1% sodium deoxycholate, 1 mM EDTA), and resuspended in 100 μL of DNA Elution buffer (50 mM sodium bicarbonate pH 8.0, 1% SDS). Sample was incubated at room temperature for 10 minutes with rotation followed by 3 minutes at 37˚C shaking, and then placed on a magnet. The supernatant was then transferred in a new tube and another 100 μL of DNA Elution buffer was added to the beads for another same incubation to obtain 200 μL of sample. Reverse-crosslink was performed by addition of Proteinase K (Qiagen, ON, Canada) and a first incubation at 55˚C for 45 minutes with shaking followed by a second incubation at 67˚C for 1.5 hours. DNA was then purified with clean up & concentrator kit (Zymo Research, Cedarlane, Canada) and quantified with QuBit (Thermo Fisher, ON, Canada). 5 μL of Streptavidin C-1 beads (Life Technologies, Thermo Fisher Scientific, ON, Canada) were washed with Tween Wash buffer (5 mM Tris-HCl pH 7.5, 0.5 mM EDTA, 1 M NaCl, 0.05% Tween-20) and resuspended in 10 μL of 2X Biotin Binding buffer (10 mM Tris-HCl pH 7.5, 1 mM EDTA, 2 M NaCl). Sample was then added for an incubation of 15 minutes at room temperature with rotation, and placed on a magnet. Supernatant was discarded and beads washed twice by adding 500 μL of Tween Wash buffer and incubating at 55˚C for 2 minutes shaking. Beads were washed in 100 μL of 1X TD buffer (2X TD buffer is 20 mM Tris-HCl pH 7.5, 10 mM magnesium chloride, 20% dimethylformamide), and then resuspended in 25 μL of 2X TD buffer, an appropriate amount of Tn5 adjusted for the DNA quantity (about 5 ng) and water to 50 μL, followed by an incubation at 55˚C for 10 minutes with interval shaking. Sample was then placed on a magnet and supernatant removed. 100 μL of 50 mM EDTA were then added to the sample for an incubation at 50˚C for 30 minutes, followed by two washes at 50˚C for 3 minutes, and two additional washes in Tween Wash buffer at 55˚C for 2 minutes. Beads were then resuspended in 10 mM Tris. Sample was placed on magnet and beads resuspended in 50 μL of a PCR master mix (25 μL of Phusion HF 2X (New England BioLabs, ON, Canada), 1 μL of each Nextera Ad1_noMX and Nextera Ad2.1 at 12.5 μM, and 23 μL of water). The following PCR program was performed: 72˚C for 5 min, 98˚C for 1 min, then cycle at 98˚C for 15 s, 63˚C for 30 s, and 72˚C for 1 min. Cycle number was set to 10 considering the amount of DNA. Library was placed on a magnet and eluted into a new tube before purification with clean up & concentrator kit (Zymo Research, Cedarlane, Canada). Library size selection (300-700 bp) was performed by cutting a 2% agarose gel and DNA purified with PureLink Quick Gel Extraction (Thermo Fisher, ON, Canada). Library quality was assessed on a BioAnalyser (Agilent, ON, Canada) using Agilent High Sensitivity DNA kit (Agilent, ON, Canada). Paired-end sequencing with 690 million reads was performed with NovaSeq 6000 S4 (UCSD IGM Genomics Facility, CA, USA). Data are available in GEO (accession number GSE154512).

H3K27ac HiChIP data processing
Paired-end reads were aligned using HiC-Pro with default settings. Contact matrix were generated with Juicer from the allValidPairs file generated by HiC-Pro [40]. Juicebox [41] was used for matrix visualisation. Most confident interactions (FDR < 0.01) were determined with FitHiChIP [42]. The matrix normalized by the vanilla coverage square root for the chromosome 6 at a 5 kb resolution was dumped using Juicer. This matrix was used to generate the virtual 4C with R considering all the interactions with a bin containing the TSS of LINC01013. H3K27ac HiChIP 1D track for VICs was generated using bamcoverage from the BAM file generated by HiC-Pro.

Chromosome conformation capture (3C)
3C experiments were performed as described previously [43]. Quantification of ligated products was performed by qPCR using the PstI cutting site located upstream the TSS of CCN2 as the anchor. For all fragments, we first normalized using a genomic GAPDH undigested control, then we used the level of ligation products between the anchor and the nearest downstream PstI site (#1) as a strong normalizer to correct for genomic random interactions as described previously [44]. Two biological replicates were conducted and qPCR were performed in triplicates to improve sensitivity. Primers sequences are provided in the S6 Table.

Cell transfection and TGF beta treatment
VICs were transfected with HiPerFect transfection reagent (Qiagen, ON, Canada). Knockdown of LINC01013 and CCN2 were performed 72 hours with a transfection of an oligonucleotide antisens (ASO) against LINC01013 (80 nM), CCN2 (40 nM) or NC5 as a negative control (IDT, IL, USA). TGF beta stimulations were performed 48 hours after 24 hours of transfection using 10 ng/mL of recombinant human TGF beta 1 (Life technologies, Thermo Fisher Scientific, ON, Canada).

Nuclear and cytoplasmic RNA purification
Nuclear and cytoplasmic quantification of LINC01013 were performed with the Cytoplasmic and Nuclear RNA Purification Kit (Norgen Biotek, ON, Canada) according to the manufacturer's instructions. U6 and HPRT were used as controls respectively for the nuclear and cytoplasmic compartments. We normalized data for each compartment with the total RNA level in cells.

Quantitative real-time PCR
RNA from mineralized and non-mineralized aortic valve tissues was purified by using the RNeasy Mini kit (Qiagen, ON, Canada). RNA from VICs was isolated with E.Z.N.A. Micro RNA kit (Omega Bio-tek, VWR, QC, Canada). One μg of RNA was reverse transcribed using the Qscript cDNA supermix from Quanta (VWR, QC, Canada). qPCR were performed with perfecta sybr supermix from Quanta on the Rotor-Gene 6000 system (Corbett Robotics Inc, CA, USA). Primers for LINC01013 were obtained from IDT (IDT, IL, USA) (F:ctaaacacaagtggtctggataga, R:tttggatcaccaaggcaag). Primers for CCN2 and COL1A1 were purchased from Qiagen (ON, Canada). The expression of the hypoxanthine-guanine phosphoribosyltransferase (HPRT) (Qiagen, ON, Canada) was used as a reference to normalize the results.

ELISA
The pro-collagen-Iα excretion by VICs was quantified from 100 μL of VICs culture medium using the Human Pro-Collagen I alpha 1 DuoSet ELISA kit (R&D systems, MN, USA) according to the manufacturer's instructions.

RIP
About 5 millions of VICs were used to perform RIP analysis. First, 20 μg of antibodies against NELF-E (Proteintech, IL, USA) or isotype IgG (Cell Signaling Technology, New England Biolabs, ON, Canada) were incubated with 50 μL of proteins G dynabeads (Life technologies, Thermo Fisher Scientific, ON, Canada) in 500 μL of IP buffer (0.5% Triton X-100, 200 mM NaCl, 10 mM Tris-HCl at pH 7.5 and 10 mM EDTA) at 4˚C overnight on a rotator. Cells were resuspended in cold lysis buffer (0.5% Triton X-100, 10 mM NaCl, 10 mM Tris-HCl pH 7.5, 10 mM EDTA, 0.5 mM PMSF, 1 mM DTT, supplemented with PIC and ProtectRNA, Sigma-Aldrich, ON, Canada) and incubated 20 minutes on ice before a centrifugation at 400 x g for 10 minutes at 4˚C. NaCl was added to a final concentration of 200 mM and lysate was precleared with 50 μL of proteins G dynabeads in IP buffer 1 hour at 4˚C on a rotator. The cleared lysate was then transferred to the tubes containing NELF-E of IgG-coated beads for 3h at 4˚C on a rotator. Samples were washed five times with 500 μL of IP buffer and then resuspended in RTK lysis buffer from E-Z Nucleic Acid (E.Z.N.A.) kit (Omega Bio-tek, VWR, QC, Canada). RNA was reverse transcribed using the Qscript cDNA supermix from Quanta (VWR, QC, Canada). qPCR were performed with perfecta sybr supermix from Quanta on the Rotor-Gene 6000 system (Corbett Robotics Inc, CA, USA).

Network analysis
Based on TissueNet [45] data, an experimentally validated dataset, protein-protein interaction (PPI) pairs were extracted in order to generate a PPI network. From data generated with the RNAs-seq in explanted mineralized AVs (TPM higher than the 10 th percentile), PPI pairs were inferred based on the TissueNet data. To generate the PPI network, the interactome of CCN2 including COL1A1, SMAD2, SMAD3 and TFGB1 was extracted. To that end, interaction pairs which included CCN2, COL1A1, SMAD2, SMAD3 and TFGB1 were isolated to form an edge list. To visualize the network, NetworkAnalyst [46] was used.

Collagen immunostaining
Human VICs were seeded on poly-L-lysine coated glass coverslips. Cells were washed once with PBS 1X and fixed in 3.7% formaldehyde for 30 minutes at room temperature. Cells were treated 20 minutes with 50 mM NH4Cl in PBS 1X and permeabilized for 10 minutes in PBS 1X containing 0.2% Triton X-100. Cells were then incubated in PBS 1X containing 5% milk for one hour at room temperature with constant agitation. Incubation with anti-COL1A1 (#E8F4L, Cell Signaling Technology, New England Biolabs, ON, Canada) was performed in PBS 1X containing 1% milk overnight at 4˚C. Cells were washed four times with PBS 1X, followed with one hour incubation with FITC-conjugated anti-rabbit secondary antibody (Molecular Probes, Thermo Fisher Scientific, ON, Canada). Secondary antibody alone was used as negative control. Slides were mounted and analysed images were acquired using a Zeiss microscope LSM800 driven by the Zen software (Objective 20X oil, 1.4 NA, Zeiss, ON, Canada). Image processing and quantification were performed with ImageJ 1.47g (NIH, USA). Integrated densities for each cell were calculated after removal of the background for the negative control.

Statistical analysis
Continuous data were expressed as mean ± SEM. Normality was tested with the Shapiro-Wilk test. For two groups, data with normal distribution were compared with Student t-test, for more than two groups with ANOVA. For data with non-normal distribution, two groups were compared with the Wilcoxon-Mann-Whitney test or with the Kruskal-Wallis test for more than two groups. A multivariate regression was used to adjust for age and sex. Categorical data were compared with Fischer's exact test. Statistical analysis were performed with JMP 14.2.0 or Prism 8.0.2.